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Spontaneous stochasticity of velocity 
in turbulence models 

Alexei A. Mailybaev* 


Abstract 

We analyze the phenomenon of spontaneous stochasticity in fluid dynamics formu¬ 
lated as the nonuniqueness of solutions resulting from viscosity at infinitesimal scales 
acting through intermediate on large scales of the flow. We study the finite-time onset 
of spontaneous stochasticity in a real version of the GOY shell model of turbulence. 
This model allows high-accuracy numerical simulations for a wide range of scales (up 
to ten orders of magnitude) and demonstrates non-chaotic dynamics, but leads to an 
infinite number of solutions in the vanishing viscosity limit after the blowup time. Thus, 
the spontaneous stochasticity phenomenon is clearly distinguished from the chaotic be¬ 
havior in turbulent flows. We provide the numerical and theoretical description of the 
system dynamics at all stages. This includes the asymptotic analysis before and after 
the blowup leading to universal (periodic and quasi-periodic) renormalized solutions, 
followed by nonunique stationary states at large times. 


1 Introduction 

One of important attributes of the developed hydrodynamic turbulence is the spontaneous 
stochasticity phenomenon. Spontaneous stochasticity is known for fields transported by tur¬ 
bulent flow [ 13 ], where non-Lipshitz inviscid velocity field leads to nonunique particle tra¬ 
jectories; see also m for a similar phenomenon for the backward in time Burgers equation. 
As the trajectories are unique at early times, when the flow is developing from large-scale 
smooth initial conditions, the hnite-time onset of spontaneous stochasticity can be expected. 
Numerical analysis of this effect is substantially limited by the extreme multiscale nature of 
the problem. The possibility of the spontaneous stochasticity in incompressible flows is closely 
related to the existence of a hnite-time blowup in the 3D Euler equations, which is the long¬ 
standing open problem of huid dynamics and applied mathematics [16] . Though there exists 
a controversy on whether the blowup was observed in numerical studies, see e.g. [smsiisoitt], 
simulations show that the blowup can be triggered by a physical boundary [ 23 ] • 

As it is common in theoretical studies, simplihed (toy) models are useful for understanding 
complex phenomena in realistic systems. One of such popular models in hydrodynamic tur¬ 
bulence is the Gledzer-Ohkitani-Yamada (GOY) shell model [T8l [3T] . modihed later into the 
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Sabra model [ 21 |. This model describes characteristic velocity variations at different scales 
considered in a geometric progression with the corresponding wave nnmbers kn = ko2"'. The 
dynamics is governed by an infinite (trnncated for nnmerical pnrposes) system of ordinary 
differential eqnations, which possesses several properties of the Navier-Stokes eqnations snch 
as scaling symmetries and conserved qnantities, and provides the chaotic behavior with the 
energy cascade between large and small scales. Thongh the GOY and Sabra models are highly 
convenient for nnmerical analysis (allowing scale resolntion over ten orders of magnitude) and 
a large amount of knowledge on their behavior is collected [3] , their dynamics is still far from 
being completely understood theoretically. 

The blowup problem for the inviscid shell models was solved [ini[2ii, showing the existence 
of a hnite-time blowup with a universal asymptotic form. A clear evidence is also given about 
the role of blowup in the turbulent energy transport to small scales dZlESlES]. In this paper, 
our goal is the numerical and theoretical description of the dynamics after the blowup in the 
vanishing viscosity limit. For this purpose, we use the real version of the (originally complex) 
GOY model. This model demonstrates the deterministic behavior before the blowup time, 
followed by spontaneous stochasticity after the blowup. Unlike its complex version, this shell 
model is not chaotic, thus, clearly distinguishing the spontaneous stochasticity from chaos in 
dynamical systems. From this perspective, the proposed model can be considered as an the 
intermediate step between the dyadic shell models ElE], whose behavior is deterministic and 
resembles the Burgers equation |29], and the chaotic complex GOY or Sabra models, whose 
inviscid limit is not yet understood. 

The contribution of this work is the description of the inviscid limit leading to an infinite 
number of solutions, triggered by vanishing viscosity through the whole range of scales. The 
specially nontrivial part of this analysis is the theory explaining the onset of spontaneous 
stochasticity, developing after the blowup. Using the renormalization technique, we show the 
existence of an asymptotic solution in the form of a universal quasi-periodic pulse moving from 
inhnitely small to large scales in finite time. The obtained results show that the nonuniqueness 
in inviscid flows with blowup can be described by considering viscosity as a small random 
variable, while keeping deterministic boundary and initial conditions. This leads to unique 
and physically relevant inviscid solutions in the form of probability measures, providing a new 
insight on the multiscale behavior of developed turbulence. 

The paper is organized as follows. We start with the model description in Section Then 
we describe viscous stationary solutions in Sectionj^and their inviscid limit in Section]^ These 
results are used in Sectionj^for defining the time-dependent inviscid solutions, which appear to 
be nonunique and depend on the specific form of inviscid limit. Section [^describes the blowup 
structure. Section [^develops the renormalization approach providing the universal asymptotic 
quasi-periodic solution after the blowup governing the onset of spontaneous stochasticity. We 
end with some conclusions. 
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2 Model and inviscid stationary solutions 


According to the Kolmogorov-Obukhov theory [211 ED], chaotic dynamics of the incompress¬ 
ible 3D Navier-Stokes equations 


du 

— + uVu = -VP 

at 


i/Au, V ■ u = 0, 


( 2 . 1 ) 


for large Reynolds numbers can be described using the concept of isotropic homogeneous tur¬ 
bulence. In this theory, the energy is transported from large (forced) to small (viscous) scales 
through a large inertial interval of intermediate scales £, where the velocity fluctuations obey 
the power-law (|5 m|) oc with the mean energy flux e from large to small scales. Such 

description, based on the dimensional analysis is known to be incomplete due to the anoma¬ 
lous corrections (not yet explained theoretically) for the exponents of this and higher velocity 
moments, see e.g. [221 [I5|. Because of the extremely large number of degrees of freedom for 
realistic Reynolds numbers, this problem is very difficult for the numerical analysis. 

In this paper, we consider the shell model of turbulence given by an inhnite system of 
differential equations 


for n = 1, 2,..., with constant boundary conditions 

Uq = const, M_i = const. 


( 2 . 2 ) 


(2.3) 


In this system, the real quantities Un model qualitatively the Fourier components of a velocity 
held u(k, t) corresponding to wavenumbers |k| = kn (one can also consider Un as a charac¬ 
teristic velocity huctuation at spatial distance £ ~ 2'K/kn), the wavenumbers kn = fco2 


are 


chosen in the form of geometric progression, and z/ is the viscosity parameter. System (2.2) 


represents the Gledzer-Ohkitani-Yamada shell model [T8l IDT] taken for purely imaginary shell 
velocities —iun] similarly, it can be deduced from the Sabra model of turbulence [ 21 ]. For 


simplicity, we set ko = 1 and consider the boundary conditions uti -|- Mq ~ 1 by an order of 


magnitude. The model (2.2) is considered as a “toy model” for 3D Navier-Stokes turbulence 


in Eq. (2.1), as it possesses the same scaling symmetry and shares some of the motion invari¬ 
ants [3|. In particular, the quadratic terms in (2.2) mimic the convective nonlinearity and the 


pressure term in Eq. (2.1). The clear advantage of the shell model is the drastic reduction in 


the number of degrees of freedom, while keeping the inhnite-dimensional nature of the how. 

Let us dehne the “energy” as P = and the “enstrophy” as D = Then the 

energy balance equation takes the form similar to the Navier-Stokes equations (see e.g. Cl) 
as 

^ = Ro - 2z/D, (2.4) 


where 


dt 


Iln knUn—lUnUn+l T 2kn-\-lUnUn-\-lUn-\-2 


(2.6) 


is the energy hux between the shells with wavenumbers kn and fcn+i- Thus, the term IIo in 


Eq. (2.4) represents the energy hux into the system generated by boundary conditions (2.3), 
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while the term 2un describes the viscous energy dissipation. A similar balance law holds for 
the helicity invariant introduced a.s H = see e.g. [3]. 

The following three transformations 


Un H- 2m„+i, z/ I-)- 2V; 

Un CUn^ U CU, t t/C] 


Un ' ^ 


il? ^n^nH-l^n+2 1? 


(2.6) 

(2.7) 

(2.8) 


are the symmetries of system ( 2 . 2 ). Also, the system is invariant under time translations 


t I—>■ t + to- Note that the shift n ^ n + 1 implies the scaling kn fcn+i = 2/c„ in Fourier 


space. Therefore, Eqs. (2.6) and (2.7) are related to the scaling of space-time. The signs an in 


Eq. (2.8) are periodic with cr„ = Un+s and can be associated with phase factors in the Fourier 


transform induced by physical space translations, see e.g. [3]. 

Stationary solutions of the inviscid model (z/ = 0) were described in [1]. The equilibrium 


conditions for Eqs. (2.2), after dividing by fcn+i, take the form 

'^n+2'^n+l 1 T g^n—l^n—2; 1,2,.... 


(2.9) 


If 7 ^ 0 for all n > 1, then Eqs. (2.9) de fine recursively the values of all shell speeds 
n > 3, for given u_i,uo,Ui,U 2 . Equations (2.9) after multiplication by 2u~\u~^i yield the 
recurrent relation for velocity ratios 




Cn — 


^n—3 


( 2 . 10 ) 


which possesses a single attracting fixed-point c„ = 1 as rz, —)■ cx). This implies that u^n+i 
is asymptotically proportional to 2~"' for each i = 1,2,3. Thus, with an increase of n, the 
stationary solutions have the period-3 asymptotic form 


u^n+i t i 1 , 2 ,3, 


( 2 . 11 ) 


with constant real quantities 01 , 02 , 03 . 


Note that the relation (2.11) implies the power law 

Un ^ l^n ) 


( 2 . 12 ) 


which can be associated with the Kolmogorov scaling law of developed turbulence mentioned 
above, as well as with the shock wave solutions in continuous models 


3 Viscous stationary solutions 


For a positive viscosity, u > 0, and constant boundary conditions (2.3), numerical simulations 


show that the shell model (2.2) possesses a fixed-point attractor. Fig. In this case, the 
equilibrium conditions become 


1 1 1 , 

Un+2Un+l T^n+l^n—1 T 2 ^ f; 2, 

4 O / 


(3.1) 
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Figure 1: Solution for the shell model for zero initial conditions, boundary values m_i = 
uq = 0.7 and viscosity z/ = 2“^^ 6 x 10“® with N = 6. Shown are the values of shell 

speeds Unit) for n = 1 , 2 ,... converging to their stationary values; the speeds decrease 
(nonmonotonously) with increasing n. 


which differ from Eq. (2.9) by an additional viscous term. For nonzero shell speeds, Eq. (3.1) 


determines the equilibrium state recursively if the four initial speeds m_i, Uq, Ui, U 2 are given. 
The speeds m_i, Uq are determined by the boundary conditions, while the speeds Ui, U 2 depend 
on the viscosity, as we will see below. 


For small viscosity z/, the viscous term in Eq. (3.1) is small and the solution is determined 


approximately by the inviscid model (2.9) leading to the properties (2.11). However, for any 


hnite z/ > 0, the viscous term becomes important for large shell numbers (small scales), when 
vkn ~ Un- Using this estimate with the scaling law (2.12), we introduce the Kolmogorov 

,- 3/4 


wavenumber as fc/c ~ z/ 


The corresponding shell number 


Uk = log2 kx = --logaZ/ 


(3.2) 


separates the region n uk of inviscid dynamics from the viscous range n > uk. Fig. 
The inviscid region contains the forcing (boundary) range at n ~ 1 and the inertial interval 
1 <C n <C Uk- 

In the viscous range, one expects fast decay of velocities with increasing shell number 


n, due to the large factor k^ of the viscous term in Eq. (2.2). In this case, the last two terms 


in Eq. (3.1) are dominant, which yields 


^n—2 ~ ^^knUn- 

Writing this expression as a linear equation 

log2 \Un\ ~ log2 \Un-l \ + log2 \Un-2\ - 2 - n - log2 V, 
one hnds the general solution in the form 

log 2 \un\ ~ ba^ + fed"' + 5 + n + log 2 z/ 


(3.3) 

(3.4) 
(3.6) 
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Figure 2 : Scaled shell speeds of the stationary attractor found numerically for the 

solution Un{t) in Fig. The viscous range position uk = 18 is shown by the vertical dashed 
line. The stationary solution is periodic in the inertial range 1 -C n -C hk- The dashed 
bold green line shows the viscous asymptotic solution kl/^^j^VnSNix) from Eq. (|4.4), which 
determines the period-3 state (4.6) in the inertial range. The inset compares numerical results 
(solid black line) with the asymptotic expression (|3.6|) (dashed blue line) for b = —0.0011; 

r 1 ' ' 

plotted are the values of log 2 | log 2 Un \ for n = 14,..., 25. 


with (T = (1 -I- -\/5)/2, cr = (1 — \/5)/2 and arbitrary real coefficients b and b. Since a > 1 and 
Id I < 1 , we hnd the asymptotic behavior as 


Ur 


32h'kr,2 




n » uk, 


(3.6) 


where the coefficient b must be negative in order to have \un\ 0 for large n (similar analysis 


was carried out in [T^ for a different shell model). Both asymptotic expressions (2.11) for 
n <C uk and (3.6) for n 3> are conhrmed numerically, see Fig. 


4 Inviscid limit for stationary solutions 

Let us denote by v^n the stationary solution corresponding to the viscosity i/, with hxed 
boundary conditions. For understanding the limit of vanishing viscosity, v —)■ -|-0, we hrst 


study the asymptotic behavior in the viscous range, n > tik, given by Eq. (3.2). 


For small viscosity, stationary solutions satisfy the period-3 asymptotic relation (2.11) in 
the inertial range 1 -C n -C uk, where the coefficients 01 , 02,03 depend on u. Expression 


( 2 . 11 ) is invariant under the scaling transformation 


Un ^ t 2 Un-\-3N i 


(4.1) 


with the shift by 3N shell numbers. The viscous range n > uk given by Eq. (3.2) can be 
shifted back to its original position, uk t uk — 3N, by changing the viscosity as 


z/ I—)■ 2 


AN 


U. 


(4.2) 
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Figure 3: Expression —)■ Vn of the limit ( 4^ shown for N = 1, 3, 5 and N ^ oo 

based on numerical simulations with the conditions of Fig. [T] 


Thus, one can expect that the combination of transformations (4.1) and (4.2) leaves the 
stationary solution approximately unchanged, i.e.. 




oN Wn] 




(4.3) 


This expression is valid for all n except for the boundary layer near n = 1, which is not shift- 


invariant. Since relation (4.3) should be satished with higher accuracy for a smaller viscosity, 
we consider the limit 


K(X) = hm = 2-"(>^+^), n, iV e Z, 


N^oo 


(4.4) 


where n is an arbitrary (positive or negative) integer and —)■ 0 as —)■ cx) for a hxed 

parameter x G M. Existence of this limit is conhrmed numerically in Fig. In Fig. |^the 
solution kn'^v^n^^ (solid line) is presented togeth er w ith its viscous range asymptotic form 

^n- 37 v^n- 3 iv(x) (dashed green line) given by Eq. (4.4). 

Combining Eqs. (2.6) and (2.7) with c = 2“^/^, we obtain the symmetry transformation 

z/ ^ 2^/=^z/, t ^ 2^/3^, (4.5) 


u. 




n-l-lj 


which leaves Eq. (2.2) invariant. Hence, expression of the limit in Eq. (4.4) represents the 


symmetry transformation (4.5) repeated 3A^ times for stationary solutions. In particular, this 


implies that 14. (y) is a stationary solution for the viscous system with v = 2“^^. The inertial 
range of I4(x) with the period-3 behavior (2.11) extends to large negative n (see Figs. |^and 
leading to the boundary condition 

yzn+i{x) Mx)k'^n+i as n ^ -cx), i = 1,2,3. (4.6) 

The coefficients Ai depend on y (see Fig. |^) and can be written using Eqs. (4.4) and (4.6) 


as 


Ai{x) = lim lim 


n—cxD V—>.oo 

m = 3(?7. -I- iV) -1- h 


i = 1,2,3. 


(4.7) 
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Figure 4: a) Functions Ai{x) determining the period-3 stationary state in the inertial interval 
depending on the parameter x for boundary conditions of Fig. b) The function a{x) 
determining the universal period-3 stationary state (4.17) in the inertial interval. Squares 
{i = 1), circles {i = 2) and diamonds {i = 3) represent the values of a = 
for n = 4 versus X ~ | obtained in numerical tests with random boundary conditions 

and random |x| < 1/2. 


The order of the limits in Eq. (4.7) is important. It describes the period-3 stationary solution 


in the inertial interval 1 ^ m -C uk, where rix = 3(x + N) —)■ oo with vanishing viscosity 


un 0. By inspection of Eq. (4.7), one can check that the relations 


Aix + t:) ^ Aix), k€Z, 


(4.8) 


are valid for any integer k. 


The obtained results show that the period-3 asymptotic form (2.11) depends on the pa¬ 


rameter X of the inviscid limit. This leads to the important conclusion that the inviscid limit 
hmy_,.+o Un^ for stationary solutions does not exist in the usual sense, while one can define an 
infinite number of inviscid stationary solutions 


Unix) = jim u 

N—^oo 




Pm = 


(4.9) 


obtained for specihc viscosity subsequences z/at —)■ 0 depending on the parameter x- These 
solutions satisfy the condition 


U3n+i(x) Mx)hnip * = 1, 2, 3, n ^ OO, 


(4.10) 


where the coefficients Aiix) are determined in Eq. (4.7) by the viscous mechanism for specihc 


boundary conditions. Since the values x and x + ^ with any integer k yield the same inviscid 
solutions, the parameter x should be considered modulo 1. 


Substituting expression (4.10) into Eq. (2.5), we hnd 


D(x) = lim n„ = 32li(x).42(x)4l3(x) 


(4.11) 
















for the energy flux in the inertial interval. This value is equal to the injected energy flux IIo 
because the viscous dissipation acts only at infinitesimal viscous scales (infinite shell numbers 
n). Since dE/dt = 0 for the stationary solution in Eq. (2.4), we find that the value D{x) 
represents the total energy dissipation rate of the inviscid solution Un{x)- This phenomenon 
of finite dissipation in the inviscid limit is known as the dissipation anomaly [T2] . 

We described the vanishing viscosity limit of stationary solutions under fixed boundary 
conditions, wh ich were taken in the numerical simulations as m_i = mq = 0.7. Using the 
symmetry (2.7) with an arbitrary coefficient c > 0, we see that Un = cUn{x) is the inviscid 


solution obtained in the vanishing viscosity limit un = c^ n = —)■ 0 for the boundary 

conditions cm_i and cuq. For this solu tion, conditions (4.10) are given by the values Ai = 
cAi{x) and the energy dissipation rate (4.11) becomes D = c^D{x)- Taking c = one 


writes z/^v = 2 Ax+^) ^^d the obtained relations can be summarized as 


^1/3 


A,; 


2-4x 


D^/^ix) 7l,(x) 2-4X’ 

The two equalities in Eq. (|4.12[) yield 


i = 1,2,3. 


(4.12) 


A- = 


Mx) ~ , loga^ , log2T>(x) 
- X+ =x +- 


L)i/3(^)’ - '12 ' 12 

By expressing y from the second equation and substituting into the first one, we obtain 


(4.13) 


Ai = X + 


loggjD 

12 


(4.14) 


where the functions ai{x) are given implicitly by the equations 


ai = 


A{x) 

des^xY 


x = x + 


log2-P(x) 

12 


(4.15) 


Numerical simulations show that the three functions ai{x) are identical up to a shift by i/3, 
i.e., (omitting tildes) 

(4.16) 


a*(x) = « ( X - 2 ) , i = 1,2,3, 


with the periodic function a{x) = a{x + 1) presented in Fig. |^. The property (4.16) can be 
understood as a result of the symmetry (4.5), which yields Ai i—)■ Aj+i and y h-)■ y — 1/3 for 


stationary solutions (4.9) 


Expressions (4.14) and (4.16) determine the period-3 asymptotic conditions (2.11) for the 


inviscid stationary solution Unix) with a given dissipation rate D as 


Usn+i dik^^n+i. 


aiix, D) = (JiD^/^a ( y - ^ ^ 


12 


i = 1, 2, 3, n —>■ oo. 


(4.17) 


where we also took into account the period-3 sign changes allowed by the symmetry (2.8). 


Conditions (4.17) are universal, i.e., valid for arbitrary boundary conditions, as we confirmed 
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numerically in Fig. for a number of tests with random values of m_i, Uq and x- Iii these 
tests, we took N = 6 corresponding to the inertial interval extended to uk = 18, and plotted 

versus y — i/3 + (log 2 F*)/12 for i = 1, 2, 3 and n = 4. 

We conclude that the viscous mechanism has the global effect on stationary solutions, 
by determining the periodic state in the inertial interval through the condition (4.17). This 


effect persists in the limit of vanishing viscosity and depends on the way this limit it taken 
governed by the parameter y. This implies the nonuniqueness of the vanishing viscosity limit, 
generating a one-parameter family of inviscid solutions, where each solution is dehned by a 
specihc viscosity sequence, upf = —)■ 0 with hxed x (mod 1). 


5 Spontaneous stochasticity in the inviscid limit 


The boundary condition ( |4.17[ ) was obtained by considering the limit (4.4), where the limiting 
expression represents the symmetry transformation (4.5) applied 3N times to a stationary 


solution. For time-dependent solutions, this transformation yields 1 1 —)■ 2‘^^t. This means that 
the time scale increases with N and, hence, the relaxation time of the time-dependent solution 


to a stationary one vanishes in the limit N —)■ oo. Therefore, the asymptotic state (4.17) must 


be attained instantaneously in the limit of vanishing viscosity z/jv = —)■ 0. 

In order to formalize this condition, let us consider a time-dependent solution Un^{t) 
dehned for viscosity z/ > 0 and hxed initial and boundary conditions. Given the parameter 
x(modl), we consider the vanishing viscosity limit as 


Unit,x) = lim u 

N^oo 


V'n] 


(«), Vk = 


(6.1) 


Then, the asymptotic condition (4.17) should hold at each time (we will see later that this 
excludes the time of blowup) as 


_1 /Q 

x)-*■ “ilx. *=1,2,3, n-t oo, 


(5.2) 


with ai{x, D) given in Eq. (4.17) and the dissipation rate D{t) depending on time. Existence 


of the limit (5.1) is conhrmed numerically and the condition (5.2) is verihed in Fig. 

Condition (5.2) implies that the vanishing viscosity limit, v —)■ -|-0, of time-dependent 
solutions does not exist in the usual sense. However, this limit can be dehned for specihc 
sequences of vanishing viscosities, z/jv —t 0. This generates an inhnite number of the inviscid 
solutions parameterized by x (modi). One can interpret such solutions in diherent ways. 
Since the viscosity v has a specihc physical value, which can be very small, the parameter x = 
— 1 log 2 z/ (mod 1) is determined uniquely. Hence, one can use this parameter for considering 
an appropriate limiting sequence z/jv = —)■ 0. Alternatively, an approximate solution 

Un{t) for small n can be found by using the viscosity z/^v with a relatively small N keeping 
the same value of x, he., one can obtain a good approximation for shell velocities at large 
scales without resolving all small scales of the model. This is analogous to the approach in 
numerical huid dynamics known as the Large Eddy Simulation, which resolves some but not 
all of the turbulence scales. 

A diherent interpretation is obtained by considering the viscosity to be small but oth¬ 
erwise undetermined (unknown). In this case, the vanishing viscosity limit may be dehned 
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Figure 5: a) Solutions for the shell u^lt) with zero initial conditions and boundary conditions 
M_i = uq = 0.7 corresponding to different values of the viscosity z/^v = The limiting 

inviscid solution depends on the parameter y after the blowup time, t > tb ^ 2.35. b) 
Dissipation rate D{t) as a function of time for the solution with x = 0 and N = 6. c) Values 
of the coefficients a* determining the solution (5.2) in the inertial interval. These coefficient 
coincide with their prediction based on the value of the dissipation rate (bold light-green 
lines) given by Eq. (4.17). 
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Figure 6 : Evolution of the shells for the solution in vanishing viscosity limit for zero 

initial conditions and boundary values m_i = Uq = 0.7. The curves show solutions for very 
small viscosities v = with y = 0, 0.04,..., 1. The inviscid solutions coincide for all 

0 < y < 1 before the blowup, t < tb, and dehne a surface for larger times, t > tb- The 
right hgure shows the amplihed region near the blowup time tb ~ 2.35, when the spontaneous 
stochasticity occurs. 


by introducing a probability distribution for the viscous parameter y as follows. Let us £x 
the boundary and initial conditions and consider y = X as a random variable uniformly 
distributed in the interval 0 < X < 1 (a random value, when chosen, is used for all times 
f > 0). Then the inviscid limit is defined as 

Unit) = lim u = (5.3) 


At each time, the solution Un\t), n = 1 , 2 ,..., represents a random variable (probability 
measure) in the U space with the norm given by the square root of energy, and the limit 
can be understood in a weak sense. Note that no choice of a special viscosity subsequence is 
necessary in the limit (5.3), where all limiting solutions (5.1) are involved through the random 
variable X. As a result, the inviscid solution Unit) is given by a singular probability measure 
supported on the one-parameter set of solutions (5.1). 

Figure shows the solution (5.3) computed numerically. One can see that the limiting 
solution is deterministic until a certain time t < tb (with the blowup time tb as described in 
the next section), and becomes stochastic for t > tb- This reveals the striking property of 
the spontaneous stochasticity of the inviscid solution Un{t) obtained in the limit of vanishing 
viscosity. 

Figur e ( thin solid line) shows the one-dimensional support of the singular probability 
measure (5.3) for the shell speeds and M 4 computed numerically at f = 3. This support 
represents a closed curve. It is remarkable that the probability measure depends on the 
viscosity mechanism. Indeed, let us consider Eq. (2.2) with the hyperviscous term —uk^Un, 
where the usual viscosity corresponds to (3 = 2. The general case with (3^2 can be studied 
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Figure 7: Shell speeds (m 3 ,M 4 ) at time t = 3 of the vanishing viscosity solu tions. Different 
viscosity mechanisms are considered given by the term —uk^Un in Eq. (2.2) for = 1.5, 2 
and 2.5. The results are obtained numerically for zero initial conditions and boundary values 
M_i = uq = 0.7. The presented curves are determined by the viscosity spanning the interval 
0 < u < uq with very small such that the viscous range starts around nx ~ 24. Inviscid 
solutions can be represented as singular probability measures (supported on a curve for the 
shells Us and U 4 ), which depend on the infinitesimal viscosity mechanism. 


in a similar way, where one should use a different scaling of viscosity depending on 13. The 
inviscid solution in the vanishing viscosity limit can be defined as a probability measure, 
similarly to Eq. (5.3). This measure was computed numerically for /3 = 1.5 and 2.5, see Fig. 
(thick solid and dashed lines). The simulations confirm that the measure is singular with 
a one-dimensional support, which is different for different f3, i.e., the inviscid limit depends 
strongly on the viscosity mechanism. 


6 Blowup 


Let us consider the Cauchy problem for the system (2.2) with arbitrary boundary conditions 
(2.3) and initial conditions at t = 0 with finite enstrophy 12(0) = ^ viscous 


case, z/ > 0, there exists a unique solution For the inviscid system, z/ = 0, the solution 

exists in a weak sense and it is unique as soon as the enstrophy is hnite. The enstrophy may 
explode (blowup) in hnite time U > 0 such that 12 —>■ 00 as f —)■ Before the blowup, 
0 < t < tb, the viscous solutions converge to the inviscid one in the vanishing viscosity limit 
z/ —)■ -|-0 |21 0 IH]. As we have shown in Section]^ the limit u +0 does not exist and the 
inviscid solution dehned for vanishing viscosity subsequences is not unique for larger times. 

The blowup structure at times t < tb was analyzed in 110] (see also [261 EU) showing that, 
at large shell numbers n, the solution has the self-similar asymptotic form 

Unit) -)■ anck~^UicCj with ^ - 4) < 0. (6.1) 

Here the scaling factor c > 0 and the signs cr^ = ±1 describe the symmetries (2.7) and (|2.8). 
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Figure 8: a) Magnified version of Fig. [^near the blowup time ti, ^ 2.35 represented by the ver¬ 
tical red line, b) Before the blowup, t < 4, the self-similar inviscid dynamics (6.1) corresponds 
to the asymptotic traveling wave ( |6.5[ ) in renormalized variables Vn{T) given by Eq. (6.4). 
Shown are the variables Vn, n = —1, 0,..., 30, at logarithmic times r = 0, Stq, . .., 21ro for 
the inviscid solution with zero initial conditions and boundary values m_i = Uq = 0.7. The 
symmetry parameters in Eq. (6.4) are c = 0.7 and (T„ = -|-1. 


The exponent y ~ 0.281 is universal (independent of boundary and initial conditions). The 
function U{^) is also universal and dehned for ^ < 0 corresponding to t < tf,, with 17(0) = 1. 
The convergence in Eq. (6.1) is pointwise for large n and arbitrary hxed i.e., for 


t = tb + ^ ^ th ^ n ^ oo. 


( 6 . 2 ) 


In our example in Fig. [^, the blowup time of the inviscid system is tb = 2.35 and one 
can recognize the self-similar behavior (6.1) for t < tb- At the blowup time, one has ^ = 0 in 
Eq. (6.1), which implies the asymptotic power law 

Unitb) —t <7nCk~'^ aS 77, —)■ CX). (6.3) 

Note that, since ?/ < 1/3, the energy flux ( 2.5[ ) tends to inhnity for large shell numbers at the 
blowup time 4, while the energy is hnite, E{tb) < oo, and the enstrophy is inhnite, Q-{tb) = oo, 
see also Fig. it- 

Let us introduce the new shell variables as 


Un = (TnCkJVn, t = tb - 2 '^/C. 


Then, following [TOl 1^ . expression (6.1) can be written in the form 


Vnir) 17 I n- 


T 

To 


where the function V{ri), rj eR, and constant tq are dehned as 

V{r]) = To = 1 - 1/^ 0.719. 


(6.4) 


(6.5) 


( 6 . 6 ) 


14 





















According to Eqs. (6.2) and (6.6), the convergence in Eq. (6.5) is understood in the limit 


T = Ton + const —)■ cx), n —>■ cx), 


(6.7) 


i.e., pointwise for a constant r] = n — t/tq. 
The function V [r]) has the limits 


lim V [i]) = 1, 

T]—^ — 00 


lim V{r]) = 0, 

TJ^OO 


( 6 . 8 ) 


where the first condition follows from the property f/(0) = 1. In the second condition, large 
Tj corresponds to the region of large shell numbers, where Un decays faster than k. 


due to 


the finite enstrophy condition and, hence, —)■ 0 in Eq. (6.4). 


Note that the limit r —)■ cx for the logarithmic time r = — log 2 [c(tb — t)] corresponds to 
t —)■ Hence, expression (6.5) describes the blowup as a traveling wave with the universal 


stationary prohle Vij]) moving from smaller to larger shell numbers with constant speed Tq”^ 
in logarithmic time r, see Fig. |^. This implies periodicity of the rescaled shell speeds n„(r) 
in Eq. (6.5), which attain the same values with the shift by one shell number after each period 
^0- 


7 Onset of spontaneous stochasticity after the blowup 

After the blowup, t > tb, the inviscid solution is not unique, as we observed in numerical 
simulations. Figs. and [^. We showed that this nonuniqueness is characterized by the 
parameter x(modl) in Eq. (5.1), determining a specihc form of the inviscid limit, = 


2 -^ix+N) g_ large times, the solutions Unit,x) converge to the stationary ones, as we 
explained in the previous sections. In this section, we focus on the behavior just after the 
blowup, i.e., in the limit f —)■ fjl", describing the onset of nonuniqueness and, thus, spontaneous 
stochasticity. 


as 


Similarly to Eqs. (6.1 )-(6.4), we introduce the functions Wn{T, y) by changing the variables 

Un = (TnCk~yWn, t = U + ‘T'" I C, (7.1) 


where c > 0 and the signs ct„ describe the symmetry transformations (2.7) and (2.8). The 


second expression introduces a new logarithmic time variable r = — log 2 [c(t — tb)] after the 
blowup, t > tb, with the right blowup limit, t —>■ corresponding to large f —)■ cx. An 
arbitrary solution Wn{f, x) can be represented as 


Wn{f, x) = ^n{n- —,X-Xc - — ) , 

\ ^0 nj 

1, 4 — 4y 

Xc = -T log 2 C, ri = ^ 2.245, 

4 1 + y 


(7.2) 


where tq = 1 — ?/ as in Eq. (6.6); the constant Xc is induced by the change of viscosity 


ujy = 2 *^x+N) symmetry transformation (2.7) with c = 2 This representation 


resembles Eq. (6.5), but takes into account the dependence on y. The reason for the choice 


of the second arguments in Wn will be clear below. 
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It is straightforward to check that the transformation 


M, 


^ fr„_ifr„2 ^u„-u t - th ^ 2^ ^(t-tb), ['(->2 


(7.3) 


is the symmetry, i.e., it relates different solntions of Eq. (2.2). It is easy to see that the 


asymptotic form (6.3) of the inviscid solntion at the blownp point {t = tb and u = 0) does not 


transforms the new variables (7.1) as 


change nnder the transformation (7.3). For inviscid solntions the symmetry (7.3) 


Wn H- Wn-u r !-)■ r + To 


with tq = 1 — y. Also, the last relation in Eq. (7.3) yields the mapping 

1 +y 


X^X + Xo, Xo = 


Ik 

Tl 


(7.4) 


(7.5) 


for the parameter of the vanishing viscosity limit = 2 —)■ 0 with ti given by Eq. (7.2). 

We see from Eqs. (7.4) and (7.5) that the argnments rji = n — f /tq and ri 2 = X ~ Xc — ^/Ti 


in Eq. (7.2) are chosen snch that the fnnctions Wn transform simply as 


Wn{r]i,r]i) H- IIA_i(771,772), 


(7.6) 


i.e., only with the change of shell nnmber. 


Recall that the asymptotic state (6.3) of the system at the blownp point tb is nniversal, i.e.. 


it is independent of initial and bonndary conditions np to the choice of the scaling parameter 
c and signs an- Hence, we can expect that similar nniversality holds after the blownp as well. 


Thns, we conjectnre (and conhrm later nnmerically) that the fnnctions in Eq. (7.2) have a 
nniversal asymptotic form for large n. This asymptotic form shonld not be affected by the 
symmetry transformation (7.3), which leaves the asymptotic state (6.3) nnchanged. Since 


this transformation changes the shell nnmber by one in the fnnctions (7.6), their nniversal 


asymptotic form, 1 ^, 1 ( 771 , 772 ) -X- 14 ^( 771 , 772 ), mnst be the independent of the shell nnmber n. 


Using Eq. (7.2), this yields the asymptotic expression of the form 


Wn{r, x) {n- —,x-Xc - — ] ■ 
To n. 


(7.7) 


Similarly to Eqs. (6.5) and (6.7), we nnderstand the limit (7.7) pointwise for large 77 —)■ cx) 


with hxed rji = n — t/tq and 112 = X~ Xc — t/ti- 

Recall that the valnes of y, which differ by an integer nnmber, correspond to the same 
inviscid solntion, see Section]^ Hence, the nniversal fnnction 11 ^( 771 , 772 ) is periodic with 
respect to the second variable as 


H"(hnh2) = W{rii,ri2 + !)• 


(7.8) 


The power law (6.3), where tb corresponds to r —)■ cx), and the relations (7.1), (7.7), (7.8) 
yield the left-side limiting valne of the fnnction IU(? 7 i, 772) with respect to the hrst argnment 
as 

lim IU( 77 i, 772 ) = 1. (7.9) 

771^—CXD 
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Figure 9: a,b) Functions Wu and Wis describing behavior after the blowup in new shell coor¬ 


dinates and logarithmic time (|7.2|). The functions are computed numerically using Eq 

,0.99 


( 2 . 2 ) 


with 35 shells for viscosities z/ = x = 0,0.01,... ,0.99. The boundary and initial 

conditions are as in Fig. c) Similar computations for the function hF 22 but for the boundary 
and initial conditions = k~'^. The functions in hgures a-c are almost identical, conhrming 
existence of the universal asymptotic form (7.7). Note that these hgures correspond to the 
same interval 3ro < f < 30ro, which yields different intervals oi t]i = n — t/tq for different n 
in Eq. (7.2). d) Renormalized shell variables Wn{f,x) computed under conditions of hgures 
a,b for X = 0 and diherent r, demonstrating quasi-periodic dynamics. 




f^n 


The limit on the other side follows from the period-3 condition (5.2), which implies that 
—)■ 0 for large n with ?/ — 1/3 < 0. This yields 


lim W{r]i,rj 2 ) = 0. 

r/i^cxD 


(7.10) 


Note that the convergence in Eq. (7.10) is rather slow due to the small absolute value of the 
exponent y — 1/3 ~ —0.0524. 


The conjectured universal asymptotic form (7.7) fully agrees with the numerical simula¬ 


tions. In Fig. |^,b we show the results of high-precision simulations carried out for Eq. (2.2) 
with the boundary conditions U-i = uq = 0.7, zero initial conditions and very small viscosities 
v = ~ 10“^^ with X = 0) 0-01) • • • > 0-99 (the viscous range starts at uk ~ 30). Shown 

are the functions hFn(hi)h 2 ) for n = lA and 18, which are determined by the corresponding 
shell speeds using Eqs. (7.1 ) and (|7.2l). These functions appear to be almost identical, conhrm¬ 


ing the asymptotic relation ( 7.7| ). Figure shows the function hF 22 (? 7 i,"^ 2 ) of the analogous 
simulation, but for the boundary and initial conditions Un = n = —1,0,1,.... These 
conditions correspond to the power law (|6.3l) at the blowup point satished exactly. Since 
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the dynamics after the blowup is quasi-periodic (together with the shift by one shell number 
after each tq), unlike the periodic dynamics before the blowup in Eq. (6.5), see Figs, [^and 
Figure shows this quasi-periodic dynamics with the logarithmic time intervals of 5tq. 
Here the viscous range starts at uk = 30, so that the formula ( 7.7| ) for the inviscid dynamics 
is appropriate for the shells n < 25. As f decreases from inhnity (the original time t increases 
in the region t > tb), the lower shell numbers get involved in the dynamics, while the higher 
shell numbers develop the period-3 tail as described in Section 

We conclude that the self-similar wave (6.5) propagating to large shell numbers (from 
large to small scales) before the blowup is followed by the wave (7.7) moving in the opposite 
direction (from small to large scales) after the blowup. The latter wave is given by a periodic 
function (7.8), which is universal (independent of initial and boundary conditions) and con¬ 
tains the information about all inviscid solutions obtained in the limit of vanishing viscosity. 
This completes the description for the onset of spontaneous stochasticity, which occurs after 
the blowup, see Fig. [^. 


8 Conclusion 

In this work, we provided the detailed numerical and theoretical description of the sponta¬ 
neous stochasticity phenomenon in the real version of the GOY shell model of turbulence. 
This deterministic model describes the energy transport from large to small scales, which is 
qualitatively similar to the one observed in the developed turbulence for the 3D Navier-Stokes 
equations. Essential reduction of the number of degrees of freedom for the shell model allows 
high-accuracy multiscale numerical analysis. We demonstrated that the limit of vanishing 
viscosity, chosen in a specihc form with a single real parameter, is not unique and leads to an 
inhnite number of different inviscid solutions. Existence of stable stationary solutions in the 
model allows understanding the inviscid limit both for stationary and time-dependent solu¬ 
tions, where inhnitesimal viscosity strongly affects the system at all scales. We showed that 
the inviscid solution, which is unique until the hnite-time blowup, becomes nonunique after the 
blowup. Emergence of the inhnite number of solutions is described using the renormalization 
technique. This approach shows that the inviscid solutions follow a universal quasi-periodic 
(in renormalized coordinates) asymptotic law triggered by the inhnitesimal viscosity, with a 
perturbation propagating from inhnitely small to large scales in hnite time. 

Considering the vanishing viscosity limit in probabilistic sense, with the viscosity as a 
small random variable, one can dehne the unique inviscid solution as a probability measure. 
We showed that this probability measure depends on the viscous mechanism: for example, it 
is diherent for a system with hyperviscosity. We stress the fact that the stochastic solution 
appears spontaneously in hnite time for an otherwise deterministic system (with no stochastic 
forcing) and for deterministic initial and boundary conditions. These results indicate the way 
for possible formalization of the spontaneous stochasticity phenomenon, which may occur in 
systems with hnite-time blowup. Note that the probability measure in our solution is singular, 
as it involves a one-dimensional set of solutions. In a forthcoming work, we give the numerical 
evidence that a similar approach leads to a continuous probability measure when applied to 
chaotic systems, i.e., the inviscid limit yields a spontaneously stochastic solution, which is 
continuously distributed in the inhnite-dimensional conhguration space. 
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